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STATE-AND- VELOCITY DEPENDENT 
FRICTIONAL SLIDING AND STRESS CORROSION DAMAGE 

CRACKING 
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Abstract. We model the progressive maturation of a heterogeneous mass to- 
wards a gravity-driven instability, characterized by the competition between 
frictional sliding and tension cracking, using axray of slider blocks on an in- 
clined basal surface, which interact via elastic-brittle springs. A realistic state- 
and rate-dependent friction law describes the block-surface interaction. The 
inner material damage occurs via stress corrosion. Three regimes, controlling 
the mass instability and its precursory behavior, are classified as a function 
of the ratio Tc/Tj of two characteristic time scales associated with internal 
damage/creep and with frictional sliding. For T^/Tf ^ 1, the whole mass un- 
dergoes a series of internal stick and slip events, associated with an initial slow 
average downward motion of the whole mass, and progressively accelerates 
until a global coherent runaway is observed. For Tc/Tj <S 1, creep/damage 
occurs sufficiently fast compared with nucleation of sliding, causing bonds to 
break, and the bottom part of the mass undergoes a fragmentation process with 
the creation of a heterogeneous population of sliding blocks. For the interme- 
diate regime Tc/Tj ~ 1, a macroscopic crack nucleates and propagates along 
the location of the largest curvature associated with the change of slope from 
the stable frictional state in the upper part to the unstable frictional sliding 
state in the lower part. The other important parameter is the Young modulus 
Y which controls the correlation length of displacements in the system. 



1. Introduction 

Gravity-driven instabilities include landslides, mountain collapses, rockfalls. ice 
mass break-off and snow avalanches. They pose a considerable risk to mountain 
communities, real-estate development, tourist activities and hydropower energy 
generation. Gravity-driven instabilities are the most widespread natural hazard 
on Earth. Most of these consist of slumping masses of soil, too small to be lethal 
but costly in term of property damage, especially when they are close to roads, 
railway lines and built-up areas. At the other extreme is the infrequent occurrence 
of the collapse of a mountainside which can release the energy equivalent of a ma- 
jor volcanic eruption or earthquake within seconds. Such collapses are among the 
most powerful hazards in nature and can wreak devastation over tens of square 
kilometers. In the US and Europe in particular, gravity-driven instabilities at all 
scales are particularly significant and cause billions of dollars and euros in damage 
each year. 
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Gravity-driven instabilities occur in a wide variety of geomechanical contexts, 
geological and structural settings, and as a response to various loading and trigger- 
ing processes. They arc often associated with other major natural disasters such 
as earthquakes, floods and volcanic eruptions. Some slope instabilities involve the 
rupture and collapse of entire mountains, with devastating consequences. One of 
the most spectacular rockslidcs that ever occurred in the world is the Kofels rock- 
slide [Erismann and Abele, 2001], in which more than 3 km'^ of rock was released in 
a short time span. This event occurred approximately 9,600 years ago. Many other 
large and giant historic and prehistoric rockslidcs have been described by Heim 
[1932], Abele [1974] and Erismann and Abele [2001]. Recent spectacular landslides 
in the European Alps include the Vaiont disaster, which resulted in about 2,000 
casualties when a 0.3 km"^ sliding mass fell into a hydro-electric power reservoir in 
1963 [Miiller, 1964]. Another example is the Val Pola slide, damming up the valley 
in 1987 [Azzoni et al., 1992]. Much can also be learned from frequent and histori- 
cally reported giant and rapid landslides on the slopes of many volcanoes [Moore et 
al., 1994; McMurtry et al., 1999; Clouard et al., 2001]. These are often associated 
with either caldera collapse [Hurlimann et al., 2000] or volcanic eruption [Voight, 
1988; Hcinrich et al, 1999]. Earthquakes are other triggering processes, which have 
to be considered. 

The present paper develops a model of the progressive maturation of a mass 
towards a gravity-driven instability, which combines basal sliding and cracking. 
Our hypothesis is that gravity-driven ruptures in natural heterogeneous material 
are characterized by a common triggering mechanism resulting from a competition 
between frictional sliding and tension cracking. Heterogeneity of material proper- 
ties and dynamical interaction seem to have a significant influence on the global 
behavior. Our main goal is thus to understand the role of the competition or inter- 
play between the two physical processes of sliding and tension cracking in the early 
stages of an instability. The run-up to the sliding instabilities can be described by 
applying a modern constitutive law of state-and-vclocity dependent friction. This 
means that solid friction is not used as a parameter but as a process evolving with 
the concentration of deformation and properties of sliding interfaces. Cracking and 
fragmentation in the mass is accounted for by using realistic laws for the progressive 
damage accumulation via stress corrosion and other thermally activated processes 
aided by stress. In one sense, the present paper can be considered as an impor- 
tant extension of [Helmstetter et al., 2004; Sornette et al., 2004], who modeled the 
potentially unstable mass by means of a single rigid slider block interacting with 
a inclined basal surface via solid friction. These authors showed that the use of a 
realistic state- and rate-dependent friction law established in the laboratory could 
provide a reasonable starting point for describing the empirical displacement and 
velocity data preceding landslides. Here, we use an array of such slider blocks on an 
inclined (and curved) basal surface, which interact via elastic-brittle springs. The 
springs, which model the inner material properties of the mass, also undergo dam- 
age, eventually leading to failure. By combining both sliding and damage processes 
in a single coherent model, this paper contributes to the literature by providing 

• a better understanding of the relative role of and interplay between the 
frictional and rupture processes, 

• a quantification of the transition from slow stable landsliding to unstable 
and fast catastrophes events. 
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• a theory of potential precursors to catastrophic slope failure, 

• a model to exploit the increasing availability of continuous monitoring by 
geodetic, GPS, SAR measurements as well as geophysical measurements to 
obtain advanced warnings, 

• a framework to assess the impact of varying climatic conditions and/or 
external forcing (rain, snow, passing seismic waves, etc.), and 

• a framework for developing strategies to mitigate the hazard resulting from 
gravity-driven instabilities. 

In Section 2, we provide a summary of previous works related to gravity-driven 
instabilities, and focus particularly on the different strategies employed to investi- 
gate stabihty. Section 3 is devoted to a detailed description of the proposed model. 
Section 4 presents an application of this numerical model to a hanging glacier and a 
comparison between numerical and experimental results. Conclusions are presented 
in Section 5. 

2. Summary of previous related works 

At the longest time scales, gravity-driven instabilities have been studied his- 
torically by the method of slope stability analysis developed in civil engineering, 
which was developed to ensure the safety of man-made structures, such as dams 
and bridges. This method usually divides a profile view of the mass of interest into 
a series of slices, and calculates the average factor of safety for each of these slices, 
using a limit equilibrium method [e.g. Hoek and Bray, 1981; Blake ct al., 2002]. 
The factor of safety is defined as the ratio of the maximum retaining force to the 
driving forces. Geomechanical data and properties are inserted in finite elements 
or discrete elements of numerical codes, which are run to predict the possible de- 
parture from static equilibrium or the distance to a failure threshold. This is the 
basis of landslide hazard maps, using safety factors greater than 1. In this modeling 
strategy, the focus is on the recognition of time-independent landslide-prone areas. 
Such an approach is similar to the practice in seismology called "time-independent 
hazard," where earthquake-prone areas are located relative to active faults, for in- 
stance, while the prediction of a single, individual earthquake is recognized to be 
much more difficult, if not unattainable. This "time-independent hazard" essen- 
tially amounts to assuming that gravity-driven instabilities result from a random 
(Poisson) process in time, and uses geomechanical modeling to constrain the future 
long-term hazard. In this model, failure appears to occur without warning, prob- 
ably because earlier movements have passed unnoticed. The approaches in terms 
of a factor of safety do not address the possible preparatory stages leading to the 
catastrophic collapse. 

In contrast, "time-dependent hazard" would accept a degree of predictability in 
the process, because the hazard associated with the potential triggering of gravity- 
driven instabilities varies with time, perhaps in association with varying external 
factors (rain, snow, volcanic). For instance, accelerated motions have been linked 
to pore pressure changes [e.g. Vangenuchten and Derijke, 1989; Van Asch et al., 
1999], leading to an instability when the gravitational pull on a slope exceeds the 
resistance at a particular subsurface level. Since pore pressure acts at the level of 
submicroscopic to macroscopic discontinuities, which themselves control the global 
friction coefficient, circulating water can hasten chemical alteration of the interface 
roughness, and pore pressure itself can force adjacent surfaces apart. Both effects 
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can produce a reduction in the friction coefficient which leads, when constant load- 
ing applies, to accelerating movements [Kuhn and Mitchell, 1993]. However, this 
explanation has not lead to a quantitative method for forecasting slope movement. 

The next level in the hierarchy of models would be "gravity-driven instability 
forecasting," which would require significantly better understanding to enable the 
prediction of some of the features of an impending catastrophe, usually on the basis 
of the observation of precursory signals. Practical difficulties include identifying 
and measuring reliable, unambiguous precursors, and the acceptance of an inherent 
proportion of missed events or false alarms. Let us also mention the approach 
developed by [Iverson, 1985; 1986a; 1986b; 1993; Iverson and Reid, 1992] based on 
coupled pressure-dependent plastic yield and nonlinear viscous flow deformations 
(see also [Iverson and Reid, 1992]). According to the theory of Iverson et al., the 
plastic yield is a generalization of the Coulomb criterion to three-dimensional stress 
states in the spirit of Rice [1975], and the effect of pore- water pressures is accounted 
for by the use of the Terzaghi "effective" stresses. This theory embodies as special 
cases models of creeping, slumping, sliding and flowing types of slope deformation. 
The underlying equations remain fundamentally based on continuum mechanics, 
which differentiates it from our approach which focuses on the dominating role of 
sliding interfaces/bands and their rheology controlled by discrete asperities, joints 
and blocks. 

There exists a common phenomenology prior to rupture, which is observed in 
many different natural materials and geophysical set-ups. In particular, there are 
well-documented cases of precursory signals, showing accelerating slip-over time 
scales of weeks to decades. We mentioned before the Vaiont landslide on the Mt 
Toe slope in the Dolomite region in the Italian Alps, which was the catastrophic 
culmination of an accelerated slope velocity (see for instance [Helmstettcr et al., 
2004; Sornette et al., 2004] for the data and a recent analysis in terms of acceler- 
ated precursory sliding velocity). More precisely, gravity-driven instabilities often 
exhibit a finite time singularity characterized by a power law behavior of some vari- 
able (such as displacement or deformation, velocity of the instable mass or acoustic 
emissions) as a function of time. This behavior is quite common in a wide class 
of nonlinear process, including gravity-driven instability such as rockfall [Amitrano 
et al., 2005], landshdes [Helmstetter et al., 2004; Sornette et al., 2004], break-off of 
ice chunks from hanging glaciers [Flotron, 1977; Rothlisberger, 1981; Luthi, 2003; 
Pralong and Funk, 2005, Pralong et al., 2005], but also for earthquakes [Bute and 
Varncs, 1993; Bowman et al., 1998; Jaumc and Sykes, 1999; Sammis and Sornette, 
2002], and volcanic eruption [Voight, 1988]. Analogous behaviors have also been 
associated with financial crashes, with the dynamics of the economy and of popu- 
lation [Johansen and Sornette, 2001; Ide and Sornette 2002]. In other cases which 
exhibited accelerated deformations , such as La Clapicre mountain [Helmstetter et 
al., 2004] in the southern Alps near St-Etienne-de-Tinee, France, the sliding can 
undergo a transition towards a slower more stable regime. While only a few such 
cases have been monitored in the past, modern monitoring techniques are bound 
to provide a wealth of new quantitative observations based on GPS and SAR (syn- 
thetic aperture radar) technology to map the surface velocity field [Mantovani et al. , 
1996; Parise, 2001; Fruneau et al., 1996; Malet et al., 2002; Bcrardino et al., 2003; 
Coe et al., 2003; Colcsanti et al., 2003; Mora et al., 2003; Wasowski and Singhroy, 
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2003] and seismic monitoring of slide quake activity [Gomberg et al., 1995; Xu et 
al., 1996; Rousseau, 1999; Caplan-Auerbach et al., 2001]. 

Other studies proposed that rates of gravity-driven mass movements are con- 
trolled by microscopic slow cracking and, when a major failure plane is developed, 
the abrupt decrease in shear resistance may provide a sufficiently large force im- 
balance to trigger a catastrophic slope rupture [Kilburn and Petley, 2003]. Such a 
mechanism, with a proper law of input of new cracks, may reproduce the accelera- 
tion preceding of the collapse that occurred at Vaiont, Mt Toe, Italy [Kilburn and 
Petley, 2003]. 

An alternative modeling strategy consists in viewing the accelerating displace- 
ment of the slope prior to the collapse as the final stage of the tertiary creep 
preceding failure [Saito and Uezawa, 1961; Saito, 1965; 1969; Voight and Kennedy, 
1979; Voight, 1988; 1989; 1990; Fukuzono, 2000; Kuhn and Mitchell, 1993]. Fur- 
ther progress in exploring the relevance of this mechanism requires a reasonable 
knowledge of the geology of the sliding surfaces, their stress-strain history, the 
mode of failure and the time-dependent shear strength along the surface of failure 
[Bhandari, 1988]. Unfortunately, this range of information is not usually available. 
This mechanism is therefore used mainly as a justification for the establishment of 
empirical criteria of impending landslide instability. 

Observations of landslides have been quantified by what amounts to a scaling law 
relating the slip acceleration dS /dt to the slip velocity S according to the following 
equation [Voight, 1990; Voight and Kennedy, 1979; Petley et al., 2002] 

(2.1) d6/dt = AS"" , 

where A and a are empirical constants. For a > 1, this relationship predicts a 
divergence of the sliding velocity in finite time at some critical time tf, which is 
determined by the parameters A and a and the initial value of the velocity. The 
mathematical divergence is of course not to be taken literally: it signals a bifurcation 
from accelerated creep to complete slope instability for which inertia is no longer 
negligible and other processes such as liquefaction (not necessarily involving fluids) 
come into play. Several cases have been quantified ex-post with this law, usually 
using a = 2, by plotting the time tf — t to failure as a function of the inverse of 
the creep velocity [Bhandari, 1988]. In fact, the solution of (|2.1[) for a > 1 can be 
expressed as 

(2.2) tf-t^ 1/6^ . 

These fits suggest that it might be possible to forecast impending landslides by 
recording accelerated precursory slope displacements. As a matter of fact, for 
Mont Toe, Italy, Vajont landslide, Voight [1988] mentioned that a prediction of the 
failure date could have been made more than 10 days before the actual failure by 
using (|2.2[) with a = 2. The physically based friction model [Helmstetter et al., 
2004; Sornette et al., 2004] which avoids such a priori assumption has confirmed 
this claim. 

Voight [1988; 1989] proposed that relation (|2.ip . which generalizes damage me- 
chanics laws [Rabotnov, 1980; Gluzman and Sornette, 2001], can be used with 
other variables (including strain and/or seismic energy release) for a large variety 
of materials and loading conditions [Sammis and Sornette, 2002]. Equation (|2.ip 
seems to apply as well to diverse types of landslides occurring in rock, soil and ice, 
including first-time and reactivated slides [Voight, 1988]. It may be seen as a special 
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Figure 1. Simplified synopsis of the range of time scales, the as- 
sociated shding velocities and the major processes involved in the 
evolution of slopes (see for example [Brunsden, 1999] for a review). 
Our model is applied to the intermediate range from years to tens 
of minutes before the catastrophic collapse, as indicated in bold- 
face. This graph was prepared with the help of G. Ouillon. 



case of a general expression for failure [Voight, 1988; 1989; Voight and Cornelius, 
1991]. Recently, such time-to- failure laws have been interpreted as resulting from 
cooperative critical phenomena and have been applied to the prediction of failure of 
heterogeneous composite materials [Anifrani et al., 1995] and to precursory increase 
of seismic activity prior to main shocks [Sornette and Sammis, 1995; Jaume and 
Sykes, 1999; Sammis and Sornette. 2002]. See also [Sornette, 2002] for extensions 
to other fields. 

The law (|2.ip may apply over time scales from years to tens of minutes before 
the run-out, as summarized in Figure [TJ Our model provides a general conceptual 
and operational framework to rationalize the law (|2.ip . One of the strengths of 
our approach is to be able to tackle in the same conceptual framework both stable 
slope evolutions as well as unstable regimes culminating in a runout. Our approach 
applies to any landslides, irrespective of their final catastrophic evolution in the 
runout regime, because we focus on the quasi-static regime in which the slope can 
still be analyzed as a system of interacting sliding blocks. While we model the 
progressive damage and rupture between sliding blocks, we do not describe the fi- 
nal stage of the dynamic runout during which most large landslide disintegrate in 
rapid motion in a collection of blocks as they move down the slope (as an example 
of this, see [Campbell 1989; 1990; Cleary and Campbell, 1993; Campbell et al., 
1995]). This final dynamic regime includes additional physical mechanisms that 
we do not account for: (i) heat-generated pore pressure inside a rapidly deforming 
shear band [Voight and Faust, 1982a, b], (ii) thermo-poro-mechanical softening of 
the soil [Vardoulakis. 2002a, b] and (ii) shear-induced pore dilation which leads to 
unstability when sediments are dilated to their critical-state porosity (see for in- 
stance [Moore and Iverson. 2002]). These models predict exponential accelerations 
only in the last tens of seconds when the sliding velocity is sufficiently large. Finally, 
the runout itself is characterized by further mechanisms associated with velocities 
reaching such large values (meter/s to tens of meters/s) that inertial effects are no 
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longer negligible (see Figure [1]) . These two times scales are not addressed in our 
approach, which focuses on the stability and possible destabilizations of slopes over 
time scales from years to tens of minutes. 

3. Geometry and physical inputs of the model 

3.1. Main elements of the model. To tackle the challenge of modeling the in- 
terplay between the two physical processes of sliding and tension cracking in the 
preparatory stages of a gravity-driven instability, we proposed a discrete set-up 
made of block masses linked to each other by elastic springs. The blocks slid on a 
basal inclined surface and interacted via direct elastic coupling to the neighboring 
blocks. 

A toy model embodying the principles of competition between subsurface slip and 
the possible growth of cracks decoupling adjacent blocks was developed by Andersen 
et al. [1997] and Leung and Andersen [1997]. Many other investigators have studied 
spring-block models, but most were derived from the Burridgc-Knopoff prototype 
having elastic coupling with an upper rigid plate, as a simple representation of the 
tectonic driving which lead to earthquakes. Note that such elastic coupling to an 
upper plane entails the wrong large scale elastic limit for gravity-driven masses. 
In contrast, the blocks in the Andersen et al. [1997] 's model are only connected 
to nearest neighbors and are in frictional contact with a substrate. Meakin and 
collaborators [Meakin, 1991] have also studied a similar model in the context of 
drying patterns. 

Using blocks as representative discrete interacting elements provides a reasonable 
starting point for describing the processes occurring before the complete destabi- 
lization of the slope. Actually, a model with just one block was used by [Heim, 
1932; Korner, 1976; Eisbacher, 1979; Davis et al., 1990]. In these works, the fric- 
tion coefficient between the rigid block and surface is usually imposed as a constant, 
or just either slip- or velocity-dependent (but not slip- and velocity-dependent). A 
constant solid friction coefficient (Mohr-Coulomb law) is often taken to simulate 
mass over bedrock sliding. A slip-dependent friction coefficient model is taken to 
simulate the yield-plastic behavior of a brittle material beyond the maximum of its 
strain-stress characteristics. For rock avalanches, Eisbacher [1979] suggested that 
the evolution from a static to a dynamic friction coefficient is induced by the emer- 
gence of a basal gauge. Studies using a velocity-dependent friction coefficient have 
focused mostly on the establishment of empirical relationships between shear stress 
T and block velocity v, such as ~ exp[ar] [Davis et al., 1990] or v ^ r^/^ [Korner, 
1976], however without a definite understanding of the possible mechanism. MuUer 
and Martel [2000] also point out that near the surface, principal stresses are either 
normal or parallel to the local slope; thus shear stress on planes parallel to the 
slope must be small. This implies that slip must initiate on pre-existing planes of 
weakness. 

Our present model improves on these models and on the multi-block model of 
Andersen et al. [1997] and Leung and Andersen [1997] in two ways. First, we use 
a state-and-velocity weakening friction law instead of a constant (or just state- or 
velocity- weakening) solid friction coefficient. Second, rather than a static threshold 
for the spring failures, we model the progressive damage accumulation via stress 
corrosion and other thermally activated processes aided by stress. Both improve- 
ments make the numerical simulations significantly more involved but present the 
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Figure 2. Illustration of the model constisting of spring-blocks 
resting on an inclined slope. The blocks lie on an inclined curved 
surface and gravity is the driving force. Only a small subset of the 
spring-block system is shown here. 

advantage of embodying rather well the known empirical phenomenology of sliding 
and damage processes. Adding the state and velocity-dependent friction law and 
time-dependent damage processes allows us to model rather faithfully the interplay 
between sliding, cracking between blocks and the overall self-organizing of the sys- 
tem of blocks. This will allow us to investigate the following questions: what are 
the conditions under which cracking (disconnection between blocks) can stabilize 
sliding? What is the effect of heterogeneity along the slope? Can we construct a 
simple set of blocks and their interactions that can reproduce the complex history 
of some slopes, such as the case studies presented in Section 

The geometry of the system of blocks interacting via springs and with a basal 
surface is depicted in Figure [51 The model includes the following characteristics: 

(1) frictional sliding on the ground or between layers, 

(2) heterogeneity of basal properties, 

(3) possible tension rupture by accumulation of damage, 

(4) dynamical interactions of damage or cracks along the sliding layer, 

(5) geometry and boundary conditions, 

(6) interplay between frictional sliding and cracking. 

We now turn to the specification of the two key ingredients, the friction and 
damage laws, that are applied to blocks and bonds respectively. 

3.2. Friction law between the discrete blocks and the basal surface . 

3.2.1. State- and rate- dependent solid friction. State- and velocity-dependent fric- 
tion laws have been established on the basis of numerous laboratory experiments 
(see, for instance, [Scholz, 2002; 1998; Marone, 1998; Gomberg et al., 2000] for 
reviews). Friction laws have been investigated in the laboratory using mainly 
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strain-controlled tests in which, for instance, motion of a slider block is driven 
by controlling its velocity. This strain-controlled set-up is thought to be relevant to 
faults which are generally assumed to be driven by far-field forces transmitted via 
elastic and inelastic deformations. In contrast, gravity-driven masses move chiefly 
in response to changes in the mix of internal body forces imposed by gravity and 
pore-pressure fields. In this context, stress-controlled experiments are more rele- 
vant. Dieterich [1994] has shown that his strain-controlled friction law does apply to 
stress-controlled interactions between ruptures. Many other works have used these 
friction law in stress-controlled situations (see e.g. [Lapusta et al., 2000; Ben-Zion 
and Rice, 1997]). Experiments have also confirmed that the state- and velocity- 
dependent friction laws established with strain-controlled experiment apply equally 
well (with a slight modification with regularization) to stress-controlled situations 
[Prakash and Clifton, 1993; Prakash, 1998; Ben-Zion, 2001; Berger, 2002]. 

Analogies between landslide faults and tectonic faults [Gomberg et al., 1995] and 
the relevance of a solid friction criterion for pre-existing slopes with weak layers 
or geologic discontinuities, with bedding planes, landslide slip-surfaces, faults, or 
joints, or for existing landslides, suggest that the improved description of the friction 
coefficient that we use here be viewed as an a clement essential to the understanding 
of slope evolutions. To our knowledge, Helmstetter et al. [2004] and Sornette et al. 
[2004] were the first to explore quantitatively the analogy between sliding rupture 
and earthquakes and to use the physics of state- and velocity-dependent friction 
to model gravity-driven instabilities and their precursory phases (see [Brabb, 1991; 
Cruden, 1997; Guzetti et al., 1999; Howell et al., 1999; Al-Homoud and Tahtamoni, 
2001] for a compilation of cases). Earlier, Chau [1995, 1999] also used the state- and 
rate-dependent friction law (|3.ip given below to model gravity-driven instabilities. 
He transformed the problem into a formal nonlinear stability analysis of the type 
found in mathematical theory of dynamical systems, but no physical application 
was described. In particular, Chau's analysis missed the existence of the finite-time 
singular behavior found also empirically, as discussed above and recalled below. 

The key concept underlying the state- and rate-dependent friction law is that 
the friction coefhcient fi between a block and the inclined basal surface supporting 
it evolves continuously as a function of cumulative slip 6, slip velocity 6 and time, 
as shown by experiments carried over a large range of time scales, under static or 
sliding conditions over a large range of velocities, for a wealth of different materials. 
The version of the rate/state- variable constitutive law, currently most accepted as 
being in reasonable agreement with experimental data on solid friction, is known 
as the Dieterich- Ruina law [Dietrich, 1994]: 

(3.1) ^i{S,0)^^io + Aln^ + BlIi■^ , 

do t'o 

where the state variable 9 is usually interpreted as the proportional to the surface 
of contact between asperities of the two surfaces. The constant is the friction 
coefficient for a sliding velocity i5o and a state variable Oq. A and B are coefficients 
depending on material properties. Expression (|3.ip as written is not appropriate 
for very low or vanishing values of sliding velocities and a simple regularization near 
5 = can be performed, that is motivated by Arrhenius thermal activation of creep 
at asperity contacts [Ben-Zion, 2003]. The time evolution of the state variable 9 is 
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(3.3) 



described by 

where Dc is a characteristic shp distance, usuaUy interpreted as the typical size of 
asperities. We note that (j3.2p can be rewritten 

For smaU or zero sliding velocity, 9 grows linearly with time. Using the above- 
mentioned regularization [Ben-Zion, 2003], this gives a logarithmic strengthening 
of /X. For a steady-state non-zero velocity, we have [Scholz, 1998] iJ. = fiQ + {A — 
B)\n\S/So], where /io = Mo + B\n-^=^. Thus, the derivative of the steady-state 
friction coefficient with respect to the logarithm of the reduced slip velocity is 
A — B. If A > B, this derivative is positive: friction increases with slip velocity 
and the system is stable as more resistance occurs which tends to react against the 
increasing velocity. In contrast, for A < B, friction exhibits the phenomenon of 
velocity-weakening and is unstable. 

The primary parameter that determines stability, A — B, is a material prop- 
erty. For instance, for granite, A ~ B is negative at low temperatures (unstable) 
and becomes positive (stable) for temperatures above about 300°C. In general, for 
low-porosity crystalline rocks, the transition from negative to positive A — B cor- 
responds to a change from elastic-brittle deformation to crystal plasticity in the 
micro-mechanics of friction [Scholz, 1998]. For the application to gravity-driven in- 
stabilities, we should in addition consider that sliding surfaces are not only contacts 
of bare rock surfaces: they are usually lined with wear detritus, called cataclasite 
or gouge in the case of faults or joints. The shearing of such granular material 
involves an additional hardening mechanism (involving dilatancy), which tends to 
make A — B more positive. For such materials, ^ — S is positive when the material 
is poorly consolidated, but decreases at elevated pressure and temperature as the 
material becomes lithified. See also Section 2.4 of Scholz's book [Scholz, 2002] and 
the discussion below. The friction law (|3.ip with (|3.2p accounts for the fundamen- 
tal properties of a broad range of surfaces in contact, namely that they strengthen 
(age) logarithmically when aging at rest, and tend to weaken (rejuvenate) when 
sliding. 

3.2.2. Frictional dynamics. Let us first consider a single block and the base in 
which it is encased. The block represents the discrete element of the slope which 
may be potentially unstable. The two main parameters controlling the stability of 
the block arc the angle cj) between the surface on which the block stands and the 
horizontal and the solid friction coefficient fi. The block exerts stresses that are 
normal (cr) as well as tangential (t) to this surface of contact. The angle (j> controls 
the ratio of the shear over normal stress: tan0 = r/a. In a first step, we assume 
for simplicity that the usual solid friction law r = fj.a holds for all times. Four 
possible regimes arc found, that depend on the ratio B/A of two parameters A and 
B of the rate and state friction law and on the initial frictional state of the sliding 
surfaces characterized by a reduced parameter 

(3.4) X, = {S 0o)'/('-^/^) ^ , fovA^B , 
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where 

(3.5) '^li^^ 

depends on the material properties but not on the initial conditions, and 9q repre- 
sents the initial value of the state variable. 

Helmstctter et al. [2004] have shown that two regimes among them can account 
for an acceleration of the displacement. For B/A > 1 (velocity weakening) and 
Xi < 1^ the slider block exhibits an unstable acceleration leading to a finite-time 
singularity of the displacement and of the velocity v ^ \/{tf ~t), thus rationalizing 
Voighfs empirical law discussed in Section 2. An acceleration of the displacement 
can also be reproduced in the velocity strengthening regime, for B/A < 1 and 
Xi > 1. In this case, the acceleration of the displacement evolves toward a stable 
sliding with a constant sliding velocity. The two others cases {B/A < 1 and Xi < 1, 
and B/A > 1 and Xi > 1) give a deceleration of the displacement. 

The case B = A is oi special interest because it retrieves the main qualitative 
features of the two classes, and also because, empirically, A is very close to i?. It 
may be natural to assume that A = B to remove the need for one parameter and 
ensure more robust results. In the sequel, we use this approximation A ^ B. In 
this case, expression p.2p renders 



j/1 

(3.6) =i-seo. 

at 

If S* 00 > 1 and is a constant, 9 decays linearly and reaches in finite time. This 
retrieves the finite-time singularity, with the slip velocity diverging as l/{tf — t) 
corresponding to a logarithmic singularity of the cumulative slip. If S Oq < 1 and 
is a constant, 6 increases linearly with time. As a consequence, the slip velocity 
decays as 5 ^ l/< at large times and the cumulative slip grows asymptotically 
logarithmically as Int. This corresponds to a standard plastic hardening behavior. 

Appendix A calculates the critical time t/ (for S9q > 1) as a function of the total 
shear force T = || X^^bond — 7weighta;|| and the normal force N = A^woight exerted 
on a given block, where i^bond is the force exerted by a neighboring spring bond, 
-^weight and Twcight arc the normal and tangential forces due to the weight of the 
block. These forces enter into the value of S defined by (|3.5p via r and a. From the 



definition of a solid friction coefficient {fi = -^), the critical time t/ signaling the 
transition from a subcritical sliding to the dynamical inertial sliding is, for /i > /iq, 

^'■'^ = exp(4-) - 1 ' 

while <j — > oo for fi < ^q. 

3.2.3. General algorithm of the frictional processes. The simulation of the frictional 
process for each given block proceeds as follows: 

(1) A given configuration of blocks and spring tensions determines the value of 
T = II X]-^bond — 7wcighta;|| and N = A^weight for each block, and therefore 
of their solid friction coefficient fj, equal to the ratio T/N. 

(2) Knowing fi for a given block and with the other material parameters l^o 
and A for that block, the time t f for the transition to the dynamical sliding 
regime for that block is calculated with expression p.7p . The formula (|3.7p 
for tf gives, therefore, the waiting time until the next slide of that block. 
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(3) When the block undergoes a transition into the dynamical sliding regime 
at time i/, its subsequent dynamics is taken to obey Newton's law. 

(4) Following other investigators [Farkas et al., 2005; Persson, 1993], we assume 
that, when the dynamical sliding phase starts, the friction coefficient de- 
creases abruptly from its value attained at time to a smaller dynamical 
value ^dyn independent of the sliding velocity. Generally, the ratio between 
the kinetic and static friction coefficients is bounded as 1/2 < ^''»"°™"= < i. 
The choice /idyn = f corresponds to the harmonic mean of the lower 
and upper bounds and is taken as a reasonable representative of the dis- 
persed set of values reported empirically. For the sake of simplicity, we 
impose this value for the dynamical friction coefficient in the dynamical 
regime. 

(5) The dynamical slide of the block goes on as long as the velocity of the block 
remains positive. When its velocity first reaches zero, we assume that the 
block stops. To account for the heterogeneity and roughness of the sliding 
surface, we assume that the state variable is reset to a new random value 
after the dynamical sliding stops. This random value is taken to reflect the 
characteristics of the new asperities constituting the fresh surface of contact. 
For the sake of simplicity, we use the same parameters do and I?c, which 
we assume identical over all blocks and constant as a function of time. 

(6) After a dynamical slide, the forces exerted by the springs that connected 
the block to its neighbors are updated, as is the new gravitational force 
(if the basal surface has a curvature), the new value of /i is obtained, the 
time counter for frictional creep is reset to zero, and a new process of slow 
frictional creep develops over the new waiting time tf, that is, in general, 
different from the previous one. 

3.3. Creep and damage process leading to tensional rupture in bonds . 

We now turn to the description of the progressive damage of the bonds between 
the blocks, which represents the possible formation of cracks and the occurrence of 
fragmentation in the mass body. 

3.3.1. Stress corrosion and creep law . It is generally known that any material sub- 
jected to a stress, constant or not, undergoes time-dependent deformation, known 
as creep. In creep, the stress is less than the mechanical strength of the mate- 
rial, but the material eventually reaches a critical state at some time tc at which a 
global catastrophic rupture occurs. By waiting a sufficiently long time, the cumu- 
lative damage building up during the deformation under stress may finally end up 
in a catastrophic rupture. Creep is all the more important, the greater the applied 
stress and the higher the temperature. The time needed to reach rupture under 
creep is controlled by the tensional versus the compressive nature of the stress and 
its magnitude, by the temperature and the microstructure of the material. Creep 
rupture phenomena have long been investigated through direct experiments [Liu 
and Ross, 1996; Guarino ct al., 2002; Lockner, 1998], and described on the basis 
of different models [Miguel et al., 2002, Ciliberto et al., 2001; Kun et al, 2003; Hi- 
dalgo et al, 2002; Main, 2000; Politi et al., 2002; Pradhan and Chakrabarti, 2003; 
Turcotte et al., 2003; Shchcrbakov and Turcotte, 2003; Vujosevic and Krajcinovic, 
1997; Saichev and Sornctte, 2005]. Many investigations focused on homogeneous 
materials such as metals [Ishikawa et al., 2002] and ceramics [Goretta et al., 2001; 
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Morita and Hiraga, 2002] and numerous recent studies are concerned with hetero- 
geneous materials such as composites and rocks [Liu and Ross, 1996; Guarino et al., 
2002; Lockner, 1998]. Visco-elastic creep and creep-rupture behaviors arc among 
the critical properties needed to assess long-term survival of material structures. 

A given bond linking two neighboring blocks is modeled as an elastic medium 
subject to creep and damage, represented by a non-linear visco-elastic rheology. 
Following Nechad et al. [2005] , the bond is assumed to be equivalent to an Eyring 
dashpot in parallel with a linear spring of stiffness E. Its deformation e is governed 
by the Eyring dashpot dynamics 

de 

(3.8) — = iirsinh(/3sdashpot) 

where the stress si in the dashpot is given by 

s 

(3.9) Sdashpot = -; 7:rrT ^ j 

1 - P(e) 

Here, s is the total stress applied to the bond and P{e) is the damage accumulated 
within the bond during its history leading to a cumulative deformation e. P(e) 
can be equivalcntly interpreted as the fraction of representative elements within 
the bond which have broken, so that the applied stress s is supported by the 
fraction 1 — P{e) of unbroken elements. The Eyring rheology (|3.8p consists, at 
the microscopic level, in adapting to the material grain rheology the theory of 
reaction rates describing processes activated by crossing potential barriers. K is a 
material property and /3 is proportional to the inverse of an effective temperature 
of activation that does not need to be the thermodynamic temperature [Ciliberto 
et al, 2001; Saichev and Sornette, 2005]. 

Following Nechad et al. [2005], we postulate the following dependence of the 
damage P(e) on the deformation e: 



where eoi,eo2 and ^ are constitutive properties of the bond material. In Nechad 
et al. [2005], expression p.lOp is derived from a mean field approximation which 
consists in assuming that the applied load is shared equally between all representa- 
tive elements (RE) in the bond: each surviving RE is subjected to the same stress 
equal to the total applied force divided by the number of surviving RE. This mean 
field approximation has been shown to be a good approximation of the elastic load 
sharing for sufficiently heterogeneous materials [Roux and Hild, 2002; Reurings and 
Alava, 2004]. The exponent ^ is a key parameter, quantifying the degree of hetero- 
geneity in the distribution of strengths of the representative elements. The smaller 
^ is, the slower P{e) goes to 1 as the deformation e increases, due to the existence 
of very strong elements still being able to support the stress. For large ^, P(e) 
moves rapidly to the value 1 corresponding to total rupture as soon as e becomes 
larger than eoi — 602- 

If the stress remains constant, it is possible to calculate the dynamics of the 
deformation e{t) following from p.8l3.9|) and ((3l0)) . Nechad et al. [2005] found 
two regimes depending on the value of s with respect to a reference value s* , defined 
as the minimum stress such that there exists at least one deformation e for which 
Sdashpot given by (|3.9|) vanishes: (i) for s < s*, the bond deformation e converges to 



(3.10) 
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a finite asymptotic value at long times and the bond does not break; (ii) for s > s* , 
the deformation e diverges in a finite time tc and the corresponding finite-time 
singularity in the damage-creep process corresponds to the rupture of the bond. 
Appendix B gives a simple formula (j5.24[) for determining tc for the bond rupture, 
when the bond is subjected to a constant stress s. 

3.3.2. Integration of stress history. The analysis of the previous subsection is tractable 
if the stress remains constant or its variation with time is simple. In our case, a 
given spring is subjected to a series of stress changes associated with the different 
sliding events of the two blocks it connects. As a consequence, the stress history 
on a given bond is in general quite complicated. 

In order to determine the time at which a bond will fail by the accumulation of 
damage, we need to integrate p.8|3.9p and (j3.10|l with a variable stress, whose time 
dependence reflects the history of sliding of the two adjacent blocks. Generalizing 
to all the bonds connecting all the blocks in the system, this leads to unwieldy 
calculations. Therefore, when the stress undergoes its first change, we use the 
following rheology to take into account arbitrary stress histories: given some stress 
history cr(<'), t' > 0, a bond is assumed to break at some fixed random time, which 
is distributed according to the following cumulative distribution function 



In other words, Fo{t) is the probability that the random time of bond failure is 
less than t. This expression amounts to the consideration of a bond failure as a 
conditional Poisson process with an intensity which is a function of all the past 
stress history weighted by the stress amplification exponent p > 0. Applied to 
material failure, this law captures the physics of failure due to stress corrosion, to 
stress-assisted thermal activation and to damage (see Newman et al. [1994; 1995] 
and references therein). The parameter p characterizes the material properties of 
the damage process. In our simulations, we will use a typical value p = 2. As 
shown by formula ()3.12p below, this value means that, if the stress is doubled, the 
waiting time until rupture is divided by four. Larger values of p would capture a 
larger sensitivity to changes in stress in this history-dependent approach. 

Consider now an element which has a given critical rupture time ti(si) given by 
equation (|5.24|) in Appendix B associated with an initial stress level si which has 
remained constant from the inception of the process till present. Let us assume 
that the stress changes from si to S2 at the time t < ti due to a change in the 
position of the blocks linked by the bond under discussion. As a consequence, the 
critical rupture time is changed from ti(si) into a new value ^12(52), which takes 
into account that some damage has already occurred until time t and that the 
damage will continue with a faster (resp. slower) rate for S2 > si (resp. S2 < Si). 
Newman et al. [1994; 1995] and Saleur et al. [1996] demonstrated that the rheology 
(|3.1ip implemented into a hierarchical system of representative elements leads to 



(3.11) 




(3.12) 



tl2{s2) = t + a{ti{si) - t) 



where 



(3.13) 
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Note that formula (|3.12p with ()3.13p recovers the fact that the critical rupture time 
remains unchanged at ti(si) if the stress remains constant at si. Note also that, 
since p > 0, a < 1 which leads to ti2(s2) < ii(si) for S2 > si, as expected. 

The new rupture time <i2(s2) (and waiting time ii2(s2) — t to rupture) can 
be compared with the rupture time t2(s2) that would have been determined by a 
constant stress S2 applied from the beginning t = 0. Indeed, for constant stresses 
si or 52 applied at t = 0, the model ()3.1ip predicts the simple scaling relation 

(3.14) Mf4^f£iy. 

^ ^ tl(si) \S2j 

For S2 > si, expression (j3.12p can be written in a form that allows a better com- 
parison with (|3.14[) as follows: 



(3 15) ^12(^2) ^ ^2(52) 

ii(si) ti(si) 



^2(52) 
ti(si) 



The scaling dependence of p.l4p is different from the exponential dependence pre- 
dicted by equation (j5.24p . that derives from the model given in the previous section 

mm 

(3.16) ^[ ^1 ^ cxp (—7(59 — Si)) , for si and S2 > s* . 

ti{si) 

However, in the limit p ^ and 7^0, both expressions (|3.14p and p.l6p become 
equivalent in the leading order of their expansions in powers of (s2 — si). In our 
simulations, we will use (|3.12p as it is the most convenient to account of all the 
sliding and rupture events that may occur in the network of blocks during the 
damage processes of all bonds. This rule embodies the creep and damage processes 
described in Section r3.3. II and by equation (|3.1ip . 

3.3.3. Possibility of bond healing. In addition, we introduce the possibility of bond 
healing, to account for self-healing properties of natural material, that may result 
from sintering and chemical reactions. Two cases have to be distinguished: 

i. In the general case, a block can be disconnected from its neighbors during the 
destabilization process. This isolated block will then slide, getting in contact with 
and pushing its neighbor situated downwards. To account for this situation, we 
assume that these blocks aggregate. In our numerical procedure, this amounts to 
removing the degrees of freedom associated with this isolated block and to doubling 
the mass of the block downward. 

ii. To prevent penetration between two blocks, we enable bonds to heal, i.e., a 
broken bond is reintroduced as soon as the distance between two blocks is smaller 
than a given threshold. 

3.3.4. Algorithm of the damage process. In summary, simulation of the damage 
process leading to bond rupture between blocks proceeds as follows. 

(1) Given an initial configuration of all the blocks within the network, the 
elastic forces exerted by all bonds can be calculated from their exten- 
sion/compression. 

(2) For each bond i subjected to an initial stress so{i), we calculate the cor- 
responding critical time tc,o(*) at which it would rupture if neither of the 
two blocks connected to it moved in the meantime. For those bonds where 
so{i) < s* defined in Equation (|5.2ip . tc,o{i) is infinite. 
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(3) However, due to the subcritical frictional processes described in Section[321 
some blocks will move when time reaches their t / given by (|3.7p . The bond 
attached to the first unstable block will have its stress modified to a new 
value in fast time during the dynamical sliding event. Note that the first 
block to slide may trigger an avalanche of block slides, therefore the stress 
is modified in more than just one bond. 

(4) For all bonds whose stresses have been changed during this "avalanche," we 
calculate their new critical rupture time using Equation (j3.12p . In general, 
each bond gets a new rupture time which is different from that of other 
bonds. 

(5) Some bonds will eventually fail, modifying the balance of forces on their 
blocks and accelerating the transition to the sliding regime, after which 
the stresses in the bonds connected to the same blocks are modified. The 
modified stresses will again lead to updates in the critical rupture times for 
those bonds, according to Expression (|3.12[) . 

(6) It is important to note that we use a separation of time scale: subcritical 
frictional creep and damage by creep preceding rupture occur in "slow time" 
compared with the "fast time" dynamics during bond rupture (which is 
assumed to occur instantaneously when the critical time tc of that bond 
is reached) and with the sliding dynamics. In other words, we freeze the 
subcritical frictional evolution as well as the damage process when blocks 
are sliding dynamically. 

This leads us now to describe the third dynamical regime occurring in our system 
of blocks. 



3.4. Dynamical interactions. When a block reaches its critical time tf for sliding 
as determined in Section 13.21 it suddenly enters a regime described by Newton's 
first equation of inertial acceleration: 



(3.17) 



Here, the total force ^ F exerted on a given block is the sum of the gravitational 
force, the dynamical frictional resistance discussed at the end of Section 13.21 and 
the elastic forces exerted by the four bonds attached to it. 

The signs and amplitudes of the four forces exerted by the four bonds depend on the 
relative position of the block with respect to the positions of its four neighboring 
blocks. Denoting by {xij , yi,j ) the coordinates of block the total force exerted 
by the four springs on this block reads 



(3.18) 
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where bij = 1 if the bond exists and bij = if the bond has failed previously and 
E is the elastic spring constant. In terms of symmetry, F^'j follows by switching 
X ^ y. 
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The dynamical friction force is opposite to the shding direction and is equal in 
amplitude to /idyn ^weight cos 4>, where Pwoight is the weight of the block and is 
the angle between the surface of contact of the block with the basal plane and the 
horizontal. 

To prevent numerical instabilities developing into oscillations close to the stopping 
point of sliding blocks, a viscosity has to be introduced in Equation (|3.17p . This 
"numerical" viscosity has a physical interpretation: radiated energy coming from 
vibrations due to jerky block sliding past asperities has been shown to lead to 
a pseudo- viscous behavior [Johansen and Sornette, 1999]. In order to avoid such 
oscillations, the viscous parameters have to be chosen in such a way that each block 
motion is in the overdamped regime. This leads to a viscous coefficient of : 



This viscous term has to be added to Equation (|3.17p . 

Note again that we freeze the subcritical frictional evolution as well as the dam- 
age process when blocks are sliding dynamically according to Equation (|3.17p . This 
allows us to decouple the fast dynamical regime described by Equation (j3.17[) from 
the subcritical frictional process occurring at the interface between each block and 
the basal surface described in Section 13.21 and the damage and creep deformation 
occurring in each bond described in Section [3.31 

In the general case, as more than one block can slip simultaneously, we will 
have to solve a system of dynamical equations like Equation (|3.17[) . one for each 
concomitantly sliding block. 

3.5. Numerical treatment and summary of our model. Each of the second- 
order differential equations in (|3.17p is transformed into a pair of first-order differen- 
tial equations. An iteration scheme using the fourth-order Runge-Kutta algorithm 
[Press at al., 1994] is adopted to solve each pair of equations. As explained above, 
all Equations (|3.17[) corresponding to the sliding blocks are solved simultaneously. 
The dynamic slip of one block may trigger the slip of other blocks, or the rupture of 
a bond. The slip of each block is computed until they all stop when their velocity 
vanishes, and the corresponding "avalanche" is terminated. Then, the subcritical 
frictional process is refreshed according to the rules of Section 13.21 and the damage 
due to creep occurring in each bond is cumulated to its previous history in line with 
the laws given in Section [3?3l 

Summary and global algorithm. The different steps in this model are described in 
Figure [3T5l As explained previously, two phases are distinguished: 

(i) A quasi-static phase corresponding to the nucleation of block sliding (Sec- 
tion [321) ^-iid bond rupture (Section [331) • 

(ii) A dynamical phase corresponding to the sliding phase of the blocks (Section 
I3.4p and of the failure of bonds. 

To account for the changes of the surface characteristics after blocks have slid, a 
random state parameter 9i is assigned to each stopping block. In this way, the 
heterogeneity of the basal properties can be reproduced and is sustained during the 
evolution of the system. 



(3.19) 




77 = 



mg cos(0) 
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Figure 3. Schematic flowchart of this modified spring-block model. 



4. Classification of the different regimes of gravity-driven 

instabilities 

111 a follow-up paper, we will apply the model presented in the previous section 
to case studies of the breaking-off of large ice masses of hanging glaciers, and in 
particular to the goal of understanding quantitatively the observations obtained 
for the Weisshorn hanging glacier [Faillettaz et al., 2008] and other glaciers in the 
Swiss Alps. 

Here, we focus on the general properties of the model. We present a classification 
of the possible regimes that emerge as a result of the interplay between the different 
mechanisms, whose amplitudes are controlled by the parameter values. 

4.1. Geometry and parameters. We consider a system of blocks, whose weights 
remain constant during the numerical simulation. This assumption is natural for 
modeling the nucleation phase preceding a rockslide or landslide. In the case of 
glaciers, this approximation is valid if the thickness of the glaciers is many times 
larger than the typical snow accumulation for a given year. For instance, in the case 
of the Weisshorn glacier with a thickness of about 30 meters, this assumption is 
well justified as, in addition, its northeastern face has a steep slope and is subjected 
to strong wind, which causes the snow to drift away. 

In order to obtain a realistic description of the damage and fragmentation that 
may develop in a sliding mass, we need a sufficiently large number of blocks. As 
a compromise between reasonable sampling and computing time, we use a model 
composed of 600 = 20 x 30 blocks for our numerical exploration of the different 
regimes. 

In this vein, consider a glacier surface area of approximately 150 m^, with a 
depth of 30 m. In a model composed of 20 x 30 blocks, each block corresponds to a 
discrete mesh of size 30 m thickness, 5 m length and 5 m width. This implies that 
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Table 1. Parameters used for the simulation. The mass lias a 
total number of Ux x Uy blocks. The second row gives the list 
of parameters. The third row specifies the units of the param- 
eters and the fourth row gives the numerical values used in the 
simulation for those parameters which have been fixed through- 
out. Missing entries correspond to the parameters that have been 
varied as explained in the text. 



the weight of each block should be approximately equal to: 0.7 x 10^ kg (with a 
density of 917 kg m"^). 

We study a particular geometry in which the slope is divided into two parts: 
Upper part stable (i.e., slope < critical slope), lower part unstable (slope > critical 
slope). Thus, the array is always driven towards an instability. Quantitatively, the 
slope (j) of the underlying supporting surface, on which the mass is sliding ranges 
from (/)top = 45° (upper part) to </>bottom = 50° (lower part). To account for the 
curvature of the slope, we used also a curved elevation of the supporting surface 
modeled by a portion of a parabola corresponding to a slope change from 45° to 50° 
over the size of the length of the glacier. We then chose a static friction parameter 
between tan 45° and tan 50°, so that the upper part is in a frictional stable state, 
while the lower part is bound to start sliding. As a consequence of the slope change, 
the sliding of the lower part will transmit deformation and stress to the upper part. 
Our goal is to study this stress and deformation transfer and document the different 
regimes of sliding, rupture and fragmentation that occur as a function of the other 
model parameters. 

The blocks are distributed in a regular mesh along the slope, in such a way that 
bonds are initially stress-free. 

The following table summarizes the parameters used in our simulations and the 
values that are fixed throughout our numerical exploration of the different regimes. 

We find that the dependence of the mass evolution as a function of the other 
parameters can be reduced to a set of two reduced parameters. 

(1) The first important parameter is the elastic coefficient E controlling the 
rigidity of the springs transferring stress between the blocks. 

(2) The second key parameter is the ratio Tc/Tf of two characteristic time scales 
associated with the two fundamental processes: internal damage/creep and 
frictional sliding. The first time scale Tc is associated with the creep-damage 
process occurring in the springs linking the blocks, also accounting for the 
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natural frequency of the springs and reads 

4.1 Te = T , 

p s* 

where T is the natural period of the spring-mass system equal to 27ry^^^^, 

and where C, (3 and s* are three parameters of the creep-damage law defined 

by expressions (|5.22p . (|5.12p and (|5.2ip . Note that 7 in expression (|5.23p 

reduces to f3 since we assume for simplicity that eo2 = eoi- 

The second time scale Tf is associated with the frictional process controlling 

the tendency of blocks to slide and is derived from expression (|5.1ip . leading 

to 

exp^ ^ ) 1 

Here, and A are two parameters of the friction law and (j)top and tan((/)bottom) 
are the two angles of the unevenly slopping basal plane supporting the in- 
stable sliding mass. 

Note that the first parameter E has also an impact on Tc through its effect on s*. 
4.2. The three main regimes. 

4.2.1. Coherent sliding of the whole mass: T^/Tf large . In this regime, the time 
needed for damage initiation is so great that the whole mass undergoes a scries 
of internal stick and slip events, associated with an initial slow average downward 
motion of the whole mass. This average motion of the mass progressively accelerates 
until a global coherent run-away is observed. Specifically, we observe the following 
evolution: 

(i) First, the blocks in the lower part, where the slope is larger, start sliding 
intermittently in a stick-slip fashion. 

(ii) The motion of these blocks transfers stress to the adjacent blocks, in turn 
leading to sliding events. 

(iii) This dynamic evolves upstream, leading to an increasing number of blocks 
sliding in a stick-slip fashion. 

(iv) The whole system is then undergoing stick-slip sliding, with increasing syn- 
chronization between the sliding blocks. 

(v) With increasing time, the system reaches a regime where a sufficiently large 
fraction of blocks are in the unstable lower part and the mass starts to slide 
as a whole, leading to a catastrophic coherent global slide. 

This sequence is illustrated in Figure 14.2.11 This regime is observed at smaller 
values of T^/Tf for small elastic coefficient E. It is observed for large values of E 
only for much larger ratios T^/Tf as shown in the phase diagram of Figure [T3l 

Figure O shows the time-dependence of four selected blocks, one in the upper 
region, two near the slope change and one in the lower region. After an initial phase, 
the velocity stabilizes and slightly increases as the whole mass slides downwards, 
where the slope is greater. 

The inset in Figure [6] shows the number of surviving bonds as a function of 
time in the upper and lower parts. One observes that damage accumulates in the 
lower unstable sliding region, whereas the upper part remains basically undamaged. 
Figure [6] shows the kinetic energy of the moving blocks as a function of time. This 
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t=2.97 1=42.37 




Figure 4. Illustration of the propagation of stick-slip events 
within the whole mass until it starts to slide coherently in a cat- 
astrophic runaway, occurring for large Tc/Tj (i.e. ^ 386), for any 
E. Time increases from left to right and top to bottom. For the 
sake of simplicity, block arc not drawn. The color of a bond evolves 
from black to red as the stress changes from lower than s* to close 
to rupture. 

graph gives an indication of the radiation energy emitted by cracks forming as a 
consequence of damage accumulation. In this regime, in which the whole mass 
evolves progressively towards a global coherent run-away sliding, one can observe 
a rich precursory behavior long before this unstable run-away occurs. 

4.2.2. Fragmentation: Tc/Tf small. This regime has characteristics which are in 
a sense opposite to those of the previous regime: creep/damage occurs before the 
nucleation of sliding has time to develop, so that bonds break and the bottom part 
of the mass undergoes a fragmentation process with the creation of a heterogeneous 
population of sliding blocks. The typical evolution of the unstable mass is as follows: 
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t (days) 



Figure 5. Displacements and veloc- 
ities of four selected blocks, one taken 
at the top of the slope (yt), two near 
slope change (j/mi and ym2) and one at 
the bottom (yt), in the stick-and-slip 
regime. 



Figure 6. Kinetic energy for the 
stick-and-slip regime. The total num- 
ber of surviving bonds is shown in the 
inset (solid line), as well as the num- 
ber of surviving bonds in the upper 
(dashed and dotted line) and lower 
part (dashed line) of the model. 



(i) First, the blocks in the lower part start sliding as in the previous case. 

(ii) As previously, the deformation due to these sliding events propagates up- 
stream, leading to stick-slip events. 

(iii) Due to the relatively fast damage process, springs between blocks break, 
further increasing the stress, promoting further spring rupture, and so on. 

(iv) Large fractured zones appear in the lower part, and blocks become isolated, 
(iv) The final state is characterized by detached isolated blocks. 

Four snapshots of the evolution of the mass illustrate this fragmentation process in 
Figure [71 This regime is observed for small values of Tc/Tf and the large elastic 
coefficient E. It is observed for small values of E, ranging up to much larger ratios 
Tc/Tf, as shown in the phase diagram of Figure [TS] 

Figure [5] shows the displacement and velocity of four selected blocks (same as 
Figure O. Note that the blocks situated in the upper part do not slide. As soon 
as the blocks start sliding, their velocity increases. Generally speaking, this can be 
explained by bond rupture (as Tc/Tf is small), meaning that the block is alone and 
does not interact anymore with its neighbors. The block accelerates according to 
equation (|3.17p with F=0. 

The figure [9] inset shows the time-dependence of the surviving bonds in the up- 
per and lower parts. A very substantial level of damage can be observed, both 
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Figure 7. Illustration of the fragmentation process associated 
with sliding, occurring for small Tc/Tf{i.c. ~ 1), for any E. Di- 
ameters of the gray circles are proportional to the corresponding 
block mass (that change after aggregation process). 



in the lower and upper parts, in contrast with the previous regime. The damage 
level in the lower part leads to a completely fragmented mass. Figure [9] plots the 
time-dependence of the kinetic energy of the moving blocks. One can observe an 
intermittent kinetic energy that increases prior to full fragmentation, thus provid- 
ing, as in the previous regime, an unambiguous precursory warning, albeit with a 
smaller advance time compared with the previous regime. 

The kinetic energy increases due to isolated blocks moving downwards after 
the rupture of their bonds (this increase starts at t—7, see Figure [7]). Because 
the system is already fragmented, this increase is of no use for an early warning. 
Radiated energy seems to provide a more reliable indication of early destabilization 
of the system. 
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Figure 8. Displacements and veloc- 
ities of four selected blocks, one taken 
at the top of the slope (yt), two near 
slope change (j/mi and ym2) and one at 
the bottom (yt), in the fragmentation 
regime. 



Figure 9. Kinetic energy in the frag- 
mentation regime. The total number 
of surviving bonds is shown in the in- 
set (solid line), as well as the num- 
ber of surviving bonds in the upper 
(dashed and dotted line) and lower 
part (dashed line) of the model. 



4.2.3. Slab avalanche: Tc/Tf ~ 1. When neither damage nor frictional sliding dom- 
inates, an interesting phenomenon is observed: the occurrence of a macroscopic 
crack propagating roughly along the location of the largest curvature associated 
with the change of slope from the stable frictional state in the upper part to the 
unstable frictional sliding state in the lower part. The evolution of the mass pro- 
ceeds as follows: 

(i) The first steps are similar to those observed in the two previous regimes, 
with individual blocks starting a stick-slip motion intermittently in the 
lower part. 

(ii) The internal stresses created by these stick-slip events lead to an upstream 
propagation of these events, which progressively invades all the blocks in 
the lower part of the mass lying on the largest slope. 

(iii) The damage time scale is such that the boost provided by the increased 
stress leads to a concentration of cracking at the boundary of the two slopes 
where the curvature is largest. 

(iv) This leads to a nucleation and growth of a macroscopic crack which sepa- 
rates the mass into two parts. 



INTERPLAY BETWEEN FRICTIONAL SLIDING AND DAMAGE CRACKING 



25 



t=1.12 t=8.22 




Figure 10. Illustration of slab avalanche process associated with 
sliding, occurring for intermediate T^/Tf (i.e. ~ 25). Diameters of 
the gray circles arc proportional to the corresponding block mass 
(that change after aggregation process) . 

(v) The lower part evolves in a run-away situation with synchronized acceler- 
ating sliding, as in the first regime described in Section 14.2.11 The upper 
part remains largely untouched and stable. 
This regime is illustrated in four snapshots shown in Figure [TOl It occurs in a band 
of large values of Tc/Tf ratios, for the small enough elastic coefBcient i5, as shown 
in the phase diagram of Figure 1131 

The Figure [T2] inset shows the number of surviving bonds as a function of time in 
the upper and lower part. Most of the damage occurs in the lower unstable sliding 
region. It can be observed that the fraction of surviving bonds drops suddenly after 
the nucleation time (here about 13 days) of the macroscopic crack separating the 
lower from the upper parts. FigurefT^ shows the time evolution of the kinetic energy 
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Figure 11. Displacements and veloc- 
ities of four selected blocks, one taken 
at the top of the slope (yt), two near 
slope change (j/mi and ym2) and one at 
the bottom (yt), in the slab avalanche 
transition regime. 



Figure 12. Kinetic energy, in the slab 
avalanche transition regime. The to- 
tal number of surviving bonds is shown 
in the inset (solid line), as well as the 
number of surviving bonds in the up- 
per (dashed and dotted line) and lower 
part (dashed line) of the model. 



of the moving blocks. Compared with previous regimes, there is a much shorter (if 
any) advance warning time, which is compatible with the first-order nature of the 
nucleation of a macroscopic crack. 

4.3. Phase diagram in the plane [E;Tc/Tf]. These three regimes are summa- 
rized in the phase diagram in the plane [E;Tc/Tf] shown in Figure [T3l 

The elastic coefficient E {= Y x I) controls the global deformation of the array 
of block-springs. Small E'-values lead to large block displacements as well as large 
deformation of the system. In contrast, for large i;^- values, the block displacements 
are small and the deformations are therefore small too. A large rigidity favors the 
emergence of a coherent behavior of the block motion. E can be understood as 
controlling the correlation length of displacements in the system. A larger iJ- value 
leads to larger clusters of blocks sliding simultaneously in local avalanches. The 
upstream propagation of the stick-slip motion is facilitated and develops faster for 
larger £"s. 

The "slab avalanche" transition regime is found for intermediate values of Tc/Tf, 
resulting from the competition of instabilities between basal friction and crack 
formations. This regime is sensitive to the presence of heterogeneities. They are 
modeled by the random resetting of the state variable 9i of the frictional process 
after each sliding event. 
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Figure 13. Phase diagram in the plane [E;Tc/Tf] summarizing 
the three main regimes identified. The transition between the 
"stick-slip transition to global runaway" and the "fi-agmentation" 
regimes is continuous. The intermediate regime is characterized 
by the occurrence of a macroscopic slab avalanche (closed to the 
unstable behaviour) or global instability (closed to the stable one). 
The existence of heterogeneities is found to play an important role 
in the transitional behavior. Fragmentation (o), stick-slip (x) and 
slab avalanche (-I-) regimes are indicated. 

5. Conclusion 

With the goal of developing a better understanding of the physical mechanisms 
leading to accelerated motions and to increasing seismic activity reported in the 
case of gravity-driven instabilities, we have developed a numerical model based on 
the discretization of the natural medium in terms of blocks and springs forming a 
two-dimensional network sliding on an inclined plane. Each block, which can slide, 
is connected to its four neighbors by springs, that can fail, depending on the history 
of displacements and damage. We develop physically realistic models describing the 
frictional sliding of the blocks on the supporting surface and the tensile failure of 
the springs between blocks proxying for crack opening. Frictional sliding is modeled 
with a state-and-velocity weakening friction law with threshold. Crack formation 
is modeled with a time-dependent cumulative damage law with thermal activation 
including stress corrosion. In order to reproduce cracking and dynamical effects, 
all equations of motion (including inertia) for each block are solved simultaneously. 

The numerical exploration of the model shows the existence of three regimes as a 
function of the ratio T^/Tf, which are the two characteristic time scales associated 
with the two fundamental processes (internal damage/creep and frictional sliding). 
For Tc/Tf ^ 1, the whole mass undergoes a series of internal stick and slip events, 
associated with an initial slow average downward motion of the whole mass which 
progressively accelerates until a global coherent run-away is observed. For Tc/Tf <IC 
1, creep/damage occurs faster than nucleation of sliding, bonds break and the 
bottom part of the mass undergoes a fragmentation process with the creation of a 
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heterogeneous population of sliding blocks. For the intermediate regime Tc/Tf ^ 1, 
a macroscopic crack forms and propagates along the location of the largest curvature 
associated with the slope change. In the upper part above the crack, the stable 
frictional state prevails, while in the lower part, the unstable frictional sliding state 
dominates. 

Our framework allows for the first time the combination of competing frictional 
and damage processes at the origin of mass instabilities. It clears the path to 
a better understanding of rupture mechanisms in heterogeneous media. It also 
casts a gleam of hope for better forecasting of the ultimate rupture of gravity- 
driven systems. Calibrating the relevant physical parameters to different natural 
materials (soil, rocks, ice or snow) will also help in the investigation of the slope 
stability under external forcing conditions (such as climate changes). 



This appendix complements Section 13.21 by providing details of the calculation 
of the critical time at which unstable sliding of a given block occurs. Section 
13.21 describes the sub-critical sliding process of a given block interacting via state- 
and-velocity solid friction with its inclined basal surface. When the sub-critical 
sliding velocity d6/dt diverges (we refer to the time when this occurs as the "critical 
time" tf for the frictional sliding instability), this signals a change of regime to the 
dynamical sliding regime where inertia (the block mass and its acceleration in the 
Newton's law) has to be taken into account. 

Let us calculate explicitly how the critical time t/ is obtained and define its de- 
pendence on the parameters and boundary conditions. Let us call T = \\Y^ ^bond — 
Twoighta;|| (or N = iVweight) the total shear (or normal) force exerted on a given 
block, where i^bond is the force exerted by a neighboring spring bond, and A^weight 
and T^^Qify ht are the normal and tangential forces due to the weight of the block. We 
then have 



where 4> is the angle of the basal surface supporting the block. Therefore, 



As explained in Section [3. 21 A — B is usually very small for natural material, of the 
order of A — B Ri ±0.02. For the sake of simplicity, we assume A ^ B. As recalled 
in Sect ion [3. 21 this choice is not restrictive as it recovers the two important regimes 
[Helmstetter et al., 2004]. This leads to 



Appendix A: Initiation of sliding for a single block 



(5.1) 



M = 




(5.2) 



/iQ + A In ^ — \- B In — = tan (p . 
So Oo 



(5.3) 




tan (j) — fiQ 



A 



whose solution is 



(5.4) 




From Equation ()3.2p . we have 



(5.5) 
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Then 

^ {So eo) exp (t^ii^) 



(5.6) e = eo+[i- ^ 

and using Equation (j5.4p . we obtain 

[So 0o) exp (t^E^^) 



(5.7) 



00 + (i - 4r exp ( ^""^-^"' ))t ' 

This expression exhibits the usual regimes: a finite time singularity is obtained for 
{So 9o) exp ( l^fL^ziia J ^ jj^ ^j^jg case, Expression (|5.7p can be re-written as 



(5.8) 

with 
(5.9) 



Dc So 00 exp 



tan (f)—fj,Q 
A 



- So 00 exp (i^Ii^) 



Dc0o 

Dc - So 00 exp (iiii^) 



We can simplify Expression (|5.8p by using the condition that, for /i = tan (/) ~ fio, 
we should have t/ ^ oo. But, for ^ = /io, exp ( ^°~^° ) — 1 and thus, for the 
condition tf s- oo to hold, we need 

(5.0) . 1 , 

The final expression for the critical time t f signaling the transition from a subcritical 
sliding to the dynamical incrtial sliding is, for /i > /ig, 

(5.11) tf = ^-^ , 

^ ^ ^ cxp(ii^)-l' 

while tf oo for fi < ^lo- Note that the dependence on So has disappeared due to 
the relation (|5T0)) . 

To summarize, a given configuration of blocks and spring tensions determines the 
value of T = 11 ^ -Fbond — ^woight^H and N = A^woight and therefore of ^ via (|5.ip . 
Knowing /i and given the other material parameters 0o,tJ-o and A, we determine 
the time tf for the transition to the dynamical regime for that block via Equation 
(ISTTD . 

Appendix B: Determination of the rupture time tc for a single bond 

Case of a constant applied stress up to rupture. Combining Equations jS 
([3:91) and ([STO]) in section [3?3l vields 



(5.12) K sinh \^^{e + eo2r - PEej 

for the equation governing the deformation e of a bond subjected to an applied 
stress s. 

Nechad et al. [2005] have studied both analytically and numerically the solution 
of this equation (|5.12p in the case s > s* (where s* is defined in the main text) 



30 



J. FAILLETTAZ, D. SORNETTE, AND M. FUNK 



for which the deformation e blows up to infinity in finite time tc, reflecting the 
accelerated tertiary creep regime culminating in the global rupture of the bond. 

Our purpose here is to provide a simple approximate formula for the rupture 
time tc, which will be used in numerical simulations of the global dynamics of the 
block array. 

An initial method to obtain tc uses an approximate expression interpolating 
between the time evolution of e{t) far from and close to tc found by Nechad et al. 
[2005]. We propose the expression 



(5.13) 



e{t) = Aln{B + Ct) -\n{- ) 



,tr 



where A,B, C, tc and r are constants to be determined from matching conditions 
far from and close to tc- In particular, Nechad et al. [2005] shows that for t <C tc, 
Equation (|5.12p simplifies into 



(5.14) 



with a = Ke 



de 

dt 



ct 



02, 



2eo2 exp(--/3s(f2a)M) and Kf3 {Eeo2 - ^J,s{^y), while for 



t tc, expression (|5.12|) simplifies into 

d 



(5.15) 



de 
'di 



- ln(<c - t) 



(I/P)-! 1 



tc-t 



with d = eoi(/3s)-i/^. 

We expand expression (|5.12|) for t ^ tc^nd for t tc and identify with Equations 
(j5.14p and (|5.15p respectively to obtain the system of five equations for the five 
unknown A,B, C, tc and r: 



(5.16) 



eo 



ac 



ac 
62 



A InB 



Ini- 



= A -Ini^ 



4 



In^^ 



^ ^ A 



c 



for f <C tc e{t = 0) = eo 



for t <^t 




for t « tc {tc - t)° 



While this system should allow us in principle to determine tc, its nonlinearity 
makes it unwieldy. We propose an alternative approach, which consists in first 
simplifying Equation (|5.13p for t <^tc and t Ki tc- 
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For i < <c 



(5.17) 



for t K tr 



de _ AC 
dt ^ B + Ct 



In 



(5.18) 



de 
'dt 



— = A IniB + Ct,)- 



/i tc — t 



In 



tr-t 



In this approach, it is consistent to put r equal to where K is the constant 
defined in equation (|3.8p . Identification with equations (|5.14|) and (|5.15p gives the 
foUowing system: 



(5.19) 



a 

b 
c 
d 



AC - In 1^ 

B , 
C , 

A hi(6 + etc) 



Ehminating A in (|5.19p yields an implicit equation that determines tc'- 



(5.20) 



-In^ 
r 



i/p 



ln(6 + etc) cd 



This implicit equation is found to be sensitive to the value of the initial parame- 
ters and its use is still rather computer-intensive, which makes it inconvenient for 
simulations involving large systems with many blocks. 

We have thus turned to a simpler but more approximate method, using the 
asymptotic analytical results obtained in Nechad et al. [2005] that, for s ^ s* 
where 



- 1\ A^-l 



is the minimum stress such that the material will fail eventually, tc is given by 

(5.22) tc = C exp(-7 s) , 
where 

(5.23) 7 = ^ • 

^01 

As K has the dimension of the inverse of time in Equation (j5.12p . we take the 
constant C to be proportional to 1/K: C ^ Our formula for determining tc in 
our simulations is therefore summarized by the following equations: 



(5.24) tc 



jr exp(— 7 s) if s > s* 
oo if s < s* 
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6. Appendix C: Energy analysis 

In order to analyze our results in term of icequake activity and compare numerical 
with experimental results, we estimate the energy evolution in the lattice. Here, 
four types of energy can be distinguished and evaluated: gravitational potential 
energy, kinetic energy, energy stored in the elastic bonds (analogous to an internal 
energy) and radiated energy. 

Gravitational potential energy. The gravitational potential energy is evaluated us- 
ing the classical expression: 



(6.1) Ep{t) ^ mg zuock 

all blocks 



Energy stored in bonds. For each bond linking block i and block j, the stored elastic 
energy is defined by: 



1 , . / ^ 2 

' E 

all bo7ids 



(6.2) Eb{t) = - 2^ E[b,,\\xj -x,\\-l 



where xj and Xi arc the position of two neighboring blocks and / the initial length 
of each bond. 

Kinetic energy. The total kinetic energy of the system is evaluated according to the 
formula 



1 /\\x(t + 6t)-x{ 



(6.3) Ek{t)^ 2"'( 



St 

all blocks 



with x{t) the position of the considered block at time t. 

Radiated energy. The radiated energy Er{t) between two time steps is evaluated 
using: 



(6.4) Er{t) ^ Et- Ep{t) - Eb{t) - Ek{t) 



with Et = -Ep(O) + Eb{0) + Ek{0) = Ep{0). In other words, it is the additional 
energy not found in the potential, elastic or kinetic energies. 
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